home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / ran2.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  3.5 KB  |  147 lines
  1. /* rng/ran2.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is an implementation of the algorithm used in Numerical
  25.    Recipe's ran2 generator.  It is a L'Ecuyer combined recursive
  26.    generator with a 32-element shuffle-box.
  27.  
  28.    As far as I can tell, in general the effects of adding a shuffle
  29.    box cannot be proven theoretically, so the period of this generator
  30.    is unknown. 
  31.  
  32.    The period of the underlying combined generator is O(2^60). */
  33.  
  34. static inline unsigned long int ran2_get (void *vstate);
  35. static double ran2_get_double (void *vstate);
  36. static void ran2_set (void *state, unsigned long int s);
  37.  
  38. static const long int m1 = 2147483563, a1 = 40014, q1 = 53668, r1 = 12211;
  39. static const long int m2 = 2147483399, a2 = 40692, q2 = 52774, r2 = 3791;
  40.  
  41. #define N_SHUFFLE 32
  42. #define N_DIV (1 + 2147483562/N_SHUFFLE)
  43.  
  44. typedef struct
  45.   {
  46.     unsigned long int x;
  47.     unsigned long int y;
  48.     unsigned long int n;
  49.     unsigned long int shuffle[N_SHUFFLE];
  50.   }
  51. ran2_state_t;
  52.  
  53. static inline unsigned long int
  54. ran2_get (void *vstate)
  55. {
  56.   ran2_state_t *state = (ran2_state_t *) vstate;
  57.  
  58.   const unsigned long int x = state->x;
  59.   const unsigned long int y = state->y;
  60.  
  61.   long int h1 = x / q1;
  62.   long int t1 = a1 * (x - h1 * q1) - h1 * r1;
  63.  
  64.   long int h2 = y / q2;
  65.   long int t2 = a2 * (y - h2 * q2) - h2 * r2;
  66.  
  67.   if (t1 < 0)
  68.     t1 += m1;
  69.  
  70.   if (t2 < 0)
  71.     t2 += m2;
  72.  
  73.   state->x = t1;
  74.   state->y = t2;
  75.  
  76.   {
  77.     unsigned long int j = state->n / N_DIV;
  78.     long int delta = state->shuffle[j] - t2;
  79.     if (delta < 1)
  80.       delta += m1 - 1;
  81.     state->n = delta;
  82.     state->shuffle[j] = t1;
  83.   }
  84.  
  85.   return state->n;
  86. }
  87.  
  88. static double
  89. ran2_get_double (void *vstate)
  90. {
  91.   float x_max = 1 - 1.2e-7f ; /* Numerical Recipes version of 1-FLT_EPS */
  92.  
  93.   float x = ran2_get (vstate) / 2147483563.0f ;
  94.  
  95.   if (x > x_max) 
  96.     return x_max ;
  97.   
  98.   return x ;
  99. }
  100.  
  101. static void
  102. ran2_set (void *vstate, unsigned long int s)
  103. {
  104.   ran2_state_t *state = (ran2_state_t *) vstate;
  105.   int i;
  106.  
  107.   if (s == 0)
  108.     s = 1;    /* default seed is 1 */
  109.  
  110.   state->y = s;
  111.  
  112.   for (i = 0; i < 8; i++)
  113.     {
  114.       long int h = s / q1;
  115.       long int t = a1 * (s - h * q1) - h * r1;
  116.       if (t < 0)
  117.     t += m1;
  118.       s = t;
  119.     }
  120.  
  121.   for (i = N_SHUFFLE - 1; i >= 0; i--)
  122.     {
  123.       long int h = s / q1;
  124.       long int t = a1 * (s - h * q1) - h * r1;
  125.       if (t < 0)
  126.     t += m1;
  127.       s = t;
  128.       state->shuffle[i] = s;
  129.     }
  130.  
  131.   state->x = s;
  132.   state->n = s;
  133.  
  134.   return;
  135. }
  136.  
  137. static const gsl_rng_type ran2_type =
  138. {"ran2",            /* name */
  139.  2147483562,            /* RAND_MAX */
  140.  1,                /* RAND_MIN */
  141.  sizeof (ran2_state_t),
  142.  &ran2_set,
  143.  &ran2_get,
  144.  &ran2_get_double};
  145.  
  146. const gsl_rng_type *gsl_rng_ran2 = &ran2_type;
  147.